IASCL Behavioral Analysis

Age at MRI and EEG

#Clean Data

All behavioral data prep

MRI correlation Data Prep

MRI Subjects, TSL LSL only

mri <- read.csv("/Users/jojohu/Documents/Qlab/iascl/RH_LH_GM_extraction.csv")

mri <- mri[,c(1,2,5,6)]

mri_gray_matter_wide <- cast(mri, part_id~roi, value = "total_gray_matter_volume")
mri_gray_matter_wide <- mri_gray_matter_wide[,!colnames(mri_gray_matter_wide) %in% "LeftAntTemp"]

average_cor_thick_wide <- cast(mri, part_id~roi, value = "average_cortical_thickness")
average_cor_thick_wide <- average_cor_thick_wide[,!colnames(average_cor_thick_wide) %in% "LeftAntTemp"]

iascl_beh_gray_matter <- 
merge(mri_gray_matter_wide, blast_spoli_data_wide,
      by.x = "part_id", by.y = "part_id",
      all.x = TRUE)

iascl_beh_gray_matter <-
  merge(iascl_beh_gray_matter, gjt[,c("part_id", "a_prime")],
        by.x = "part_id", by.y = "part_id",
        all.x = TRUE)

iascl_beh_gray_matter <-
iascl_beh_gray_matter[,c(1:12, 25, 36, 16, 17, 22, 23, 27, 28, 31, 32)]

iascl_beh_gray_matter <- 
merge(iascl_beh_gray_matter, drive_checklist[,c("Participant.ID", "age_at_neuro_year")],
      by.x = "part_id", by.y = "Participant.ID",
      all.x = TRUE)

iascl_beh_gray_matter_td <-
iascl_beh_gray_matter[which(iascl_beh_gray_matter$group == "TD"),] 

iascl_beh_gray_matter_asd <-
iascl_beh_gray_matter[which(iascl_beh_gray_matter$group == "ASD"),] 


iascl_beh_cor_thick <- 
merge(average_cor_thick_wide, blast_spoli_data_wide,
      by.x = "part_id", by.y = "part_id",
      all.x = TRUE)

iascl_beh_cor_thick <-
  merge(iascl_beh_cor_thick, gjt[,c("part_id", "a_prime")],
        by.x = "part_id", by.y = "part_id",
        all.x = TRUE)

iascl_beh_cor_thick <-
iascl_beh_cor_thick[,c(1:12, 25, 36, 16, 17, 22, 23, 27, 28, 31, 32)]

iascl_beh_cor_thick <- 
merge(iascl_beh_cor_thick, drive_checklist[,c("Participant.ID", "age_at_neuro_year")],
      by.x = "part_id", by.y = "Participant.ID",
      all.x = TRUE)

iascl_beh_cor_thick_td <-
iascl_beh_cor_thick[which(iascl_beh_cor_thick$group == "TD"),] 

iascl_beh_cor_thick_asd <-
iascl_beh_cor_thick[which(iascl_beh_cor_thick$group == "ASD"),] 

#MRI data participants
print(iascl_beh_cor_thick[,c("part_id", "group")])
##        part_id group
## 1  blast_c_019    TD
## 2  blast_c_059    TD
## 3  blast_c_061   ASD
## 4  blast_c_071   ASD
## 5  blast_c_095   ASD
## 6  blast_c_126    TD
## 7  blast_c_147   ASD
## 8  blast_c_165    TD
## 9  blast_c_168    TD
## 10 blast_c_206    TD
## 11 blast_c_207    TD
## 12 blast_c_210    TD
## 13 blast_c_213    TD
## 14 blast_c_215  <NA>
## 15 blast_c_216    TD
## 16 blast_c_224    TD
#Intersection between Beh and MRI participants
intersect(blast_spoli_data_wide$part_id, iascl_beh_cor_thick$part_id)
##  [1] "blast_c_019" "blast_c_059" "blast_c_061" "blast_c_071" "blast_c_095"
##  [6] "blast_c_126" "blast_c_147" "blast_c_165" "blast_c_168" "blast_c_206"
## [11] "blast_c_207" "blast_c_210" "blast_c_213" "blast_c_216" "blast_c_224"
identical(iascl_beh_cor_thick$part_id, iascl_beh_gray_matter$part_id)
## [1] TRUE

EEG + MRI participants, SSL, TSL only

Task, group two-way Ancova controlling for age

#Extract participants that have all SL task completed
iascl_all_completed <- 
  eeg_mri_beh_data[,c(1:7,11:12,17:18,20,35,22,23,26,27,37:38)]

iascl_all_completed_acc <- iascl_all_completed[,c(1, 6, 2, 8, 9)]
iascl_all_completed_slope <- iascl_all_completed[,c(1, 6, 2, 16, 17)]

iascl_all_completed_acc<- melt(iascl_all_completed_acc, 
                               id.vars = c("part_id", "group", "age_at_neuro_year"))
## Warning: attributes are not identical across measure variables; they will
## be dropped
iascl_all_completed_slope<- melt(iascl_all_completed_slope, 
                               id.vars = c("part_id", "group", "age_at_neuro_year"))
## Warning: attributes are not identical across measure variables; they will
## be dropped
colnames(iascl_all_completed_acc)[colnames(iascl_all_completed_acc) == "variable"] <- "task"
colnames(iascl_all_completed_acc)[colnames(iascl_all_completed_acc) == "value"] <- "accuracy"

colnames(iascl_all_completed_slope)[colnames(iascl_all_completed_slope) == "variable"] <- "task"
colnames(iascl_all_completed_slope)[colnames(iascl_all_completed_slope) == "value"] <- "rt_slope"
#aov(accuracy~group*task + age+Error(subjid/task), data = data)
acc_task_group_aov <- 
aov(accuracy~group*task + age_at_neuro_year + Error(part_id/task), data = iascl_all_completed_acc)
## Warning in aov(accuracy ~ group * task + age_at_neuro_year + Error(part_id/
## task), : Error() model is singular
summary(acc_task_group_aov)
## 
## Error: part_id
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## group              1 0.0321 0.03205   2.169 0.1533  
## task               1 0.0230 0.02301   1.557 0.2237  
## age_at_neuro_year  1 0.0457 0.04565   3.089 0.0911 .
## group:task         1 0.0415 0.04153   2.810 0.1061  
## Residuals         25 0.3695 0.01478                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: part_id:task
##            Df  Sum Sq  Mean Sq F value Pr(>F)
## task        1 0.00195 0.001946   0.196  0.663
## group:task  1 0.00008 0.000082   0.008  0.929
## Residuals  18 0.17873 0.009929
rt_slope_task_group_aov <- 
aov(rt_slope~group*task + age_at_neuro_year + Error(part_id/task), data = iascl_all_completed_slope)
## Warning in aov(rt_slope ~ group * task + age_at_neuro_year + Error(part_id/
## task), : Error() model is singular
summary(rt_slope_task_group_aov)
## 
## Error: part_id
##                   Df   Sum Sq   Mean Sq F value Pr(>F)  
## group              1 0.001114 0.0011136   3.116 0.0898 .
## task               1 0.000036 0.0000363   0.102 0.7526  
## age_at_neuro_year  1 0.000204 0.0002044   0.572 0.4566  
## group:task         1 0.000001 0.0000007   0.002 0.9643  
## Residuals         25 0.008936 0.0003574                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: part_id:task
##            Df   Sum Sq   Mean Sq F value Pr(>F)
## task        1 0.000748 0.0007482    1.01  0.328
## group:task  1 0.000852 0.0008522    1.15  0.298
## Residuals  18 0.013333 0.0007407

Gender Chi-square

iascl_all_completed_gender_wide <- iascl_all_completed[,c("part_id", "group", "sex")]

iascl_all_completed[,c(2:5, 8:19)] <- lapply(iascl_all_completed[,c(2:5, 8:19)], function(x) {
  if (is.factor(x))
    as.numeric(as.character(x))
  else
    x
})
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
iascl_all_completed_gender_long <- melt(iascl_all_completed_gender_wide, id.vars = c("part_id", "group"))

gender_wide <- cast(iascl_all_completed_gender_long, group~value, length)

# Chi-square on gender
# p<0.05 tells us that they are not matched for gender between group
print(gender_wide)
##   group  F  M
## 1   ASD  1 11
## 2    TD 14  8
print(chisq.test(gender_wide))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  gender_wide
## X-squared = 7.5199, df = 1, p-value = 0.006102

Gender is NOT matched between group.

Age difference between group

td_age <- 
  iascl_all_completed[which(iascl_all_completed$group == "TD"), 
                       c("part_id", "group", "age_at_neuro_year")]

td_age$age_at_neuro_month <- td_age$age_at_neuro_year *12

asd_age <- 
  iascl_all_completed[which(iascl_all_completed$group == "ASD"), 
                       c("part_id", "group", "age_at_neuro_year")]

asd_age$age_at_neuro_month <- asd_age$age_at_neuro_year *12

library("pastecs")
## 
## Attaching package: 'pastecs'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## The following objects are masked from 'package:data.table':
## 
##     first, last
td_age_descrip_stat <- stat.desc(td_age)
td_age_descrip_stat <- td_age_descrip_stat[c("nbr.val","mean", "std.dev"),"age_at_neuro_year"]
td_age_descrip_stat <- data.frame(td_age_descrip_stat)
td_age_descrip_stat <- cbind(c("n","mean", "std.dev"), td_age_descrip_stat)
print(td_age_descrip_stat)
##   c("n", "mean", "std.dev") td_age_descrip_stat
## 1                         n           22.000000
## 2                      mean            9.284136
## 3                   std.dev            2.037441
asd_age_descrip_stat <- stat.desc(asd_age)
asd_age_descrip_stat <- asd_age_descrip_stat[c("nbr.val","mean", "std.dev"),"age_at_neuro_year"]
asd_age_descrip_stat <- data.frame(asd_age_descrip_stat)
asd_age_descrip_stat <- cbind(c("n","mean", "std.dev"), asd_age_descrip_stat)
print(asd_age_descrip_stat)
##   c("n", "mean", "std.dev") asd_age_descrip_stat
## 1                         n           12.0000000
## 2                      mean            7.9305833
## 3                   std.dev            0.9043391
#Test age difference between group, p <0.05 means age is significantly different between group
print(t.test(td_age$age_at_neuro_year, asd_age$age_at_neuro_year))
## 
##  Welch Two Sample t-test
## 
## data:  td_age$age_at_neuro_year and asd_age$age_at_neuro_year
## t = 2.6708, df = 31.151, p-value = 0.01192
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.3201399 2.3869662
## sample estimates:
## mean of x mean of y 
##  9.284136  7.930583

Age IS significantly different between TD and ASD (for all in lab participants.)

Prep data for graphs, remove outliers

iascl_sl_bar_all <- iascl_all_completed[,c(1:6, 8:11, 14:17)]


#Remove outliners for accuracy data
ssl_accurarcy_outlier <-
iascl_sl_bar_all[which(iascl_sl_bar_all$accuracy_children_ssl_accuracies < 0.15), "part_id"]

tsl_accurarcy_outlier <-
iascl_sl_bar_all[which(iascl_sl_bar_all$accuracy_children_tsl_accuracies < 0.15), "part_id"]

if(length(ssl_accurarcy_outlier) != 0) {
  iascl_sl_bar_all[which(
    iascl_sl_bar_all$part_id %in% ssl_accurarcy_outlier),
  c("accuracy_children_ssl_accuracies")] <- NA
}

if(length(ssl_accurarcy_outlier) != 0) {
  iascl_all_completed[which(
    iascl_all_completed$part_id %in% ssl_accurarcy_outlier),
  c("accuracy_children_ssl_accuracies")] <- NA
}


if(length(tsl_accurarcy_outlier) != 0) {
  iascl_sl_bar_all[which(
    iascl_sl_bar_all$part_id %in% tsl_accurarcy_outlier),
  c("accuracy_children_tsl_accuracies")] <- NA
}

if(length(tsl_accurarcy_outlier) != 0) {
  iascl_all_completed[which(
    iascl_all_completed$part_id %in% tsl_accurarcy_outlier),
  c("accuracy_children_tsl_accuracies")] <- NA
}

#Remove outliners for RT and RT slope data-----------------------------------------------------------------
ssl_rt_hit_count <- read.csv("ssl_rt_hit_trial_count.csv")
tsl_rt_hit_count <- read.csv("tsl_rt_hit_trial_count.csv")

ssl_rt_hit_count <- 
  ssl_rt_hit_count[which(ssl_rt_hit_count$part_id %in% iascl_sl_bar_all$part_id),]
tsl_rt_hit_count <- 
  tsl_rt_hit_count[which(tsl_rt_hit_count$part_id %in% iascl_sl_bar_all$part_id),]

ssl_rt_outlier <- ssl_rt_hit_count[which(ssl_rt_hit_count$hit_trial_number < 6), "part_id"]
tsl_rt_outlier <- tsl_rt_hit_count[which(tsl_rt_hit_count$hit_trial_number < 6), "part_id"]


if(length(ssl_rt_outlier) != 0) {
  iascl_sl_bar_all[which(
    iascl_sl_bar_all$part_id %in% ssl_rt_outlier),
  c("rt_children_ssl_indiv_rts",
    "slope_children_ssl_indiv_rts_slope")] <- NA
}

if(length(ssl_rt_outlier) != 0) {
  iascl_all_completed[which(
    iascl_all_completed$part_id %in% ssl_rt_outlier),
  c("rt_children_ssl_indiv_rts",
    "slope_children_ssl_indiv_rts_slope")] <- NA
}

if(length(tsl_rt_outlier) != 0) {
  iascl_sl_bar_all[which(
    iascl_sl_bar_all$part_id %in% tsl_rt_outlier),
  c("rt_children_tsl_indiv_rts",
    "slope_children_tsl_indiv_rts_slope")] <- NA
}


if(length(tsl_rt_outlier) != 0) {
  iascl_all_completed[which(
    iascl_all_completed$part_id %in% tsl_rt_outlier),
  c("rt_children_tsl_indiv_rts",
    "slope_children_tsl_indiv_rts_slope")] <- NA
}

These are removed outliers:

SSL RT and Slope < 6 hits:

TSL RT and Slope < 6 hits: blast_c_209

SSL Accuracy:

TSL Accuracy:

Group SL analysis & KBITT, SCQ Descriptive Stat

T-test TD vs. ASD

td_sl <- 
  iascl_all_completed[which(iascl_all_completed$group == "TD"),]


asd_sl <- 
  iascl_all_completed[which(iascl_all_completed$group == "ASD"),]


library("pastecs")
td_kbit_scq_descrip_stat <- stat.desc(td_sl[,c(3:4,5)])
td_kbit_scq_descrip_stat <- td_kbit_scq_descrip_stat[c("nbr.val","mean", "std.dev"),
                                                     c("kbit_matrices_raw", "kbit_matrices_std", "scq_total")]
print(td_kbit_scq_descrip_stat)
##         kbit_matrices_raw kbit_matrices_std scq_total
## nbr.val         22.000000          22.00000 20.000000
## mean            31.409091         113.18182  2.150000
## std.dev          6.870635          16.01379  1.565248
asd_kbit_scq_descrip_stat <- stat.desc(asd_sl[,c(3:4,5)])
asd_kbit_scq_descrip_stat <- asd_kbit_scq_descrip_stat[c("nbr.val","mean", "std.dev"),
                                                     c("kbit_matrices_raw", "kbit_matrices_std", "scq_total")]
print(asd_kbit_scq_descrip_stat)
##         kbit_matrices_raw kbit_matrices_std scq_total
## nbr.val         12.000000          12.00000 11.000000
## mean            26.333333         106.91667 13.909091
## std.dev          8.370221          18.54948  5.048852
#p <0.05, TD and ASD are truly different
#t.test(td_sl$accuracy_children_ssl_accuracies, asd_sl$accuracy_children_ssl_accuracies)

t_test_multi_pair <- 
  function(x,y){
  test <- t.test(x,y)

  data.frame(p_value = test$p.value,
             df = test$parameter,
             t_stat = test$statistic)
  }


#Results withOUT outliers
sapply(intersect(colnames(td_sl),colnames(asd_sl))[c(8:11, 14:17)], 
                 function(x) t_test_multi_pair(td_sl[,x], asd_sl[,x]))
##         accuracy_children_ssl_accuracies accuracy_children_tsl_accuracies
## p_value 0.5655127                        0.135311                        
## df      17.27191                         24.96517                        
## t_stat  0.5859128                        1.543421                        
##         entropy_children_ssl_entropy entropy_children_tsl_entropy
## p_value 0.9006048                    0.24506                     
## df      20.72792                     18.14556                    
## t_stat  -0.1264365                   -1.201355                   
##         rt_children_ssl_indiv_rts rt_children_tsl_indiv_rts
## p_value 0.02945021                0.2593682                
## df      15.33708                  23.84279                 
## t_stat  2.400955                  1.155416                 
##         slope_children_ssl_indiv_rts_slope
## p_value 0.1201148                         
## df      16.13298                          
## t_stat  -1.641148                         
##         slope_children_tsl_indiv_rts_slope
## p_value 0.8822999                         
## df      19.88359                          
## t_stat  -0.149969
sapply(intersect(colnames(td_sl),colnames(asd_sl))[c(3:4, 5)], 
                 function(x) t_test_multi_pair(td_sl[,x], asd_sl[,x]))
##         kbit_matrices_raw kbit_matrices_std scq_total   
## p_value 0.08818368        0.3356349         1.117523e-05
## df      19.2119           20.02764          11.0689     
## t_stat  1.796338          0.9865457         -7.52821

No group difference is seen for KBIT.

SSL Mean RT is significantly different between group (TD > ASD). No outliners in the in lab sample.

SCQ significantly higher in ASD.

Ancova comparing between group for all measures

ancova_bt_group <- 
  function(x,y,z){
  test <- aov(x~y + age_at_neuro_year, z)
  summary(test)
  }


#Results with outliers
sapply(colnames(iascl_all_completed)[c(8:11, 14:17)], 
       function(x) ancova_bt_group(iascl_all_completed[,x], iascl_all_completed[,"group"], iascl_all_completed))
## $accuracy_children_ssl_accuracies
##                   Df   Sum Sq   Mean Sq F value Pr(>F)
## y                  1 0.003110 0.0031097  0.3363 0.5684
## age_at_neuro_year  1 0.005692 0.0056921  0.6156 0.4419
## Residuals         20 0.184924 0.0092462               
## 
## $accuracy_children_tsl_accuracies
##                   Df  Sum Sq  Mean Sq F value  Pr(>F)  
## y                  1 0.03767 0.037673  2.2743 0.14459  
## age_at_neuro_year  1 0.05279 0.052792  3.1870 0.08687 .
## Residuals         24 0.39755 0.016565                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $entropy_children_ssl_entropy
##                   Df   Sum Sq   Mean Sq F value Pr(>F)
## y                  1 0.000163 0.0001633  0.0131 0.9101
## age_at_neuro_year  1 0.000030 0.0000299  0.0024 0.9614
## Residuals         20 0.249840 0.0124920               
## 
## $entropy_children_tsl_entropy
##                   Df  Sum Sq  Mean Sq F value  Pr(>F)  
## y                  1 0.03307 0.033065  1.2618 0.27242  
## age_at_neuro_year  1 0.16940 0.169400  6.4644 0.01787 *
## Residuals         24 0.62892 0.026205                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $rt_children_ssl_indiv_rts
##                   Df Sum Sq Mean Sq F value  Pr(>F)  
## y                  1 162142  162142  5.9346 0.02432 *
## age_at_neuro_year  1   6211    6211  0.2273 0.63867  
## Residuals         20 546428   27321                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $rt_children_tsl_indiv_rts
##                   Df Sum Sq Mean Sq F value Pr(>F)
## y                  1  37580   37580  1.2691 0.2715
## age_at_neuro_year  1  47880   47880  1.6170 0.2162
## Residuals         23 681049   29611               
## 
## $slope_children_ssl_indiv_rts_slope
##                   Df    Sum Sq    Mean Sq F value Pr(>F)
## y                  1 0.0015102 0.00151017  2.8766 0.1054
## age_at_neuro_year  1 0.0008359 0.00083595  1.5923 0.2215
## Residuals         20 0.0104998 0.00052499               
## 
## $slope_children_tsl_indiv_rts_slope
##                   Df    Sum Sq    Mean Sq F value Pr(>F)
## y                  1 0.0000103 0.00001028  0.0225 0.8820
## age_at_neuro_year  1 0.0000105 0.00001052  0.0230 0.8807
## Residuals         23 0.0104977 0.00045642

Ancova controlling for age shows two groups are still significantly different in SSL mean RT. Also, Age had an significant effect on TSL entropy (F = 6.4644, p = 0.01787). Age on TSL accuracy (F = 3.1870, p = 0.08687)

Accuracy against chance for TD and ASD

t_test_against_chance <- 
  function(x){
  test <- t.test(x, mu=0.5, alternative= "greater")

  data.frame(p_value = test$p.value,
             df = test$parameter,
             t_stat = test$statistic)
  }

#TD without outliers
sapply(colnames(td_sl)[c(8:9)], 
       function(x) t_test_against_chance(td_sl[,x]))
##         accuracy_children_ssl_accuracies accuracy_children_tsl_accuracies
## p_value 0.0515594                        0.0111262                       
## df      13                               15                              
## t_stat  1.753091                         2.548765
#ASD without outliers
sapply(colnames(asd_sl)[c(8:9)], 
       function(x) t_test_against_chance(asd_sl[,x]))
##         accuracy_children_ssl_accuracies accuracy_children_tsl_accuracies
## p_value 0.2640249                        0.2747999                       
## df      8                                10                              
## t_stat  0.6595811                        0.6192483

TD: marginally significant above chance for SSL accuracy (t = 1.75, p = 0.051) Significantly above chance for TSL accuracy (t = 2.55, p = 0.011)

ASD: Not above chance for either.

Ignore: Accuracy bar graph without outliers

Accuracy Line Plot (no outliers)

detach(package:plyr)
library(dplyr)
library(ggplot2)

#The group_by is strange in dplyr, n keeps changing when I change mutate labels
# iascl_sl_bar_acc %>% 
#   filter(value != "")  %>% 
#   mutate(group = factor(group, labels = c("ASD", "TD")), 
#          variable = factor(variable, labels = c("Tone", "Syllable"))) %>% 
#   group_by(group, variable) %>%
#   summarise(acc_mean = mean(value, na.rm = TRUE),
#             acc_se = sd(value, na.rm = TRUE)/ sqrt(n()),
#             n = n())

#Descriptive Stats for TD Tone Accuracy
library("pastecs")
stat_desc <- 
function(df,group,task) {
  stat.desc(df[which(df$group == group & df$variable == task),])
}

td_tone <- stat_desc(iascl_sl_bar_acc, "TD", "Tone")
td_syllable <- stat_desc(iascl_sl_bar_acc, "TD", "Syllable")
asd_tone <- stat_desc(iascl_sl_bar_acc, "ASD", "Tone")
asd_syllable <- stat_desc(iascl_sl_bar_acc, "ASD", "Syllable")
 
extract_line_graph_stat <- function (x, group, task) {
y <- x[c("nbr.val","mean", "SE.mean"),"value"]
dim(y) <- c(3, 1)
y <- t(y)
colnames(y) <- c("n", "mean", "std_error")
y <- data.frame(y)
y$group <- group
y$task <- task


return(y)
}

td_tone_stat <- extract_line_graph_stat(td_tone, "TD", "Tone")
td_syllable_stat <- extract_line_graph_stat(td_syllable, "TD", "Syllable")
asd_tone_stat <- extract_line_graph_stat(asd_tone, "ASD", "Tone")
asd_syllable_stat <- extract_line_graph_stat(asd_syllable, "ASD", "Syllable")

acc_stat <- rbind(td_tone_stat, td_syllable_stat, asd_tone_stat, asd_syllable_stat)

colnames(acc_stat)[colnames(acc_stat) == "group"] <- "Group"

library(plyr)
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:reshape':
## 
##     rename, round_any
acc_stat$task <-
revalue(acc_stat$task, c(#"entropy_children_lsl_entropy"="Letter",
                                     #"entropy_children_vsl_entropy"="Image",
                                     "Syllable"="Linguistic",
                                     "Tone"="Non-linguistic"))

pd <- position_dodge(width = 0.2)

acc_stat %>%
  ggplot(aes(x = task, y = mean, group = Group)) +
  geom_line(aes(linetype = Group, color = Group), position = pd, size = 1.8) +
   geom_errorbar(aes(ymin = mean - std_error, ymax = mean + std_error), width = .1, position = pd) +
     geom_point(aes(color = Group), size = 4, position = pd,show.legend = FALSE) + 
          geom_point(size = 3, color = "white", position = pd) +
      labs(title = "Behavioral Accuracy across Groups",
         x = "SL Task",  # Change x-axis label
         y = "Mean Accuracy") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(
        plot.title = element_text(size=16, face="bold"),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text=element_text(size=12, face = "bold")
        ) +
  theme(legend.text=element_text(size=14, face="bold"),
        legend.title=element_text(size=15, face="bold")) +
   theme(
    panel.background = element_rect(fill = "white"),         # Set plot background to white
    legend.key  = element_rect(fill = "white"),              # Set legend item backgrounds to white
    axis.line.x = element_line(colour = "black", size = 1),  # Add line to x axis
    axis.line.y = element_line(colour = "black", size = 1)   # Add line to y axis
  ) +
  geom_hline(yintercept=0.5, color = "grey") +
  geom_text(aes(x= 2.3, y = 0.505, label = paste0("chance", "-level")), size=5) +
  geom_text(aes(x= 2.05, y = 0.64, label = paste0("*")), size=7) +
  geom_text(aes(x= 1.05, y = 0.64, label = paste0("†")), size=7) 

# ggsave("behavioral_acc_across_group.png",
#        width = 15, height = 15, units = "cm")

Ignore: Entropy bar graph without outliers

Entropy Line Plot (no outliers)

detach(package:plyr)
library(dplyr)
library(ggplot2)

library("pastecs")
stat_desc <- 
function(df,group,task) {
  stat.desc(df[which(df$group == group & df$variable == task),])
}

td_tone <- stat_desc(iascl_sl_bar_entropy, "TD", "Tone")
td_syllable <- stat_desc(iascl_sl_bar_entropy, "TD", "Syllable")
asd_tone <- stat_desc(iascl_sl_bar_entropy, "ASD", "Tone")
asd_syllable <- stat_desc(iascl_sl_bar_entropy, "ASD", "Syllable")
 
extract_line_graph_stat <- function (x, group, task) {
y <- x[c("nbr.val","mean", "SE.mean"),"value"]
dim(y) <- c(3, 1)
y <- t(y)
colnames(y) <- c("n", "mean", "std_error")
y <- data.frame(y)
y$group <- group
y$task <- task


return(y)
}

td_tone_stat <- extract_line_graph_stat(td_tone, "TD", "Tone")
td_syllable_stat <- extract_line_graph_stat(td_syllable, "TD", "Syllable")
asd_tone_stat <- extract_line_graph_stat(asd_tone, "ASD", "Tone")
asd_syllable_stat <- extract_line_graph_stat(asd_syllable, "ASD", "Syllable")

entropy_stat <- rbind(td_tone_stat, td_syllable_stat, asd_tone_stat, asd_syllable_stat)

pd <- position_dodge(width = 0.2)

entropy_stat %>%
  ggplot(aes(x = task, y = mean, group = group)) +
  geom_line(aes(linetype = group), position = pd) +
   geom_errorbar(aes(ymin = mean - std_error, ymax = mean + std_error), width = .1, position = pd) +
     geom_point(size = 4, position = pd) + 
   geom_point(size = 3, color = "white", position = pd) +
   guides(linetype = guide_legend("Group")) +
      labs(title = "Behavioral SL Task Entropy",
         x = "Task",  # Change x-axis label
         y = "Entropy") +
  theme(plot.title = element_text(hjust = 0.5)) + #Center title
   theme(
    panel.background = element_rect(fill = "white"),         # Set plot background to white
    legend.key  = element_rect(fill = "white"),              # Set legend item backgrounds to white
    axis.line.x = element_line(colour = "black", size = 1),  # Add line to x axis
    axis.line.y = element_line(colour = "black", size = 1)   # Add line to y axis
  ) 

Mean RT bar graph without outliers

Mean RT Line Plot (no outliers)

detach(package:plyr)
library(dplyr)
library(ggplot2)

library("pastecs")
stat_desc <- 
function(df,group,task) {
  stat.desc(df[which(df$group == group & df$variable == task),])
}

td_tone <- stat_desc(iascl_sl_bar_rt, "TD", "Tone")
td_syllable <- stat_desc(iascl_sl_bar_rt, "TD", "Syllable")
asd_tone <- stat_desc(iascl_sl_bar_rt, "ASD", "Tone")
asd_syllable <- stat_desc(iascl_sl_bar_rt, "ASD", "Syllable")
 
extract_line_graph_stat <- function (x, group, task) {
y <- x[c("nbr.val","mean", "SE.mean"),"value"]
dim(y) <- c(3, 1)
y <- t(y)
colnames(y) <- c("n", "mean", "std_error")
y <- data.frame(y)
y$group <- group
y$task <- task


return(y)
}

td_tone_stat <- extract_line_graph_stat(td_tone, "TD", "Tone")
td_syllable_stat <- extract_line_graph_stat(td_syllable, "TD", "Syllable")
asd_tone_stat <- extract_line_graph_stat(asd_tone, "ASD", "Tone")
asd_syllable_stat <- extract_line_graph_stat(asd_syllable, "ASD", "Syllable")

rt_stat <- rbind(td_tone_stat, td_syllable_stat, asd_tone_stat, asd_syllable_stat)

pd <- position_dodge(width = 0.2)

rt_stat %>%
  ggplot(aes(x = task, y = mean, group = group)) +
  geom_line(aes(linetype = group), position = pd) +
   geom_errorbar(aes(ymin = mean - std_error, ymax = mean + std_error), width = .1, position = pd) +
     geom_point(size = 4, position = pd) + 
   geom_point(size = 3, color = "white", position = pd) +
   guides(linetype = guide_legend("Group")) +
      labs(title = "Behavioral SL Task Mean Reaction Time",
         x = "Task",  # Change x-axis label
         y = "Mean Reaction Time (ms)") +
  theme(plot.title = element_text(hjust = 0.5)) + #Center title
   theme(
    panel.background = element_rect(fill = "white"),         # Set plot background to white
    legend.key  = element_rect(fill = "white"),              # Set legend item backgrounds to white
    axis.line.x = element_line(colour = "black", size = 1),  # Add line to x axis
    axis.line.y = element_line(colour = "black", size = 1)   # Add line to y axis
  ) 

Scaled RT Slope bar graph without outliers

RT Slope Line Plot (no outliers)

detach(package:plyr)
library(dplyr)
library(ggplot2)

library("pastecs")
stat_desc <- 
function(df,group,task) {
  stat.desc(df[which(df$group == group & df$variable == task),])
}

td_tone <- stat_desc(iascl_sl_bar_slope, "TD", "Tone")
td_syllable <- stat_desc(iascl_sl_bar_slope, "TD", "Syllable")
asd_tone <- stat_desc(iascl_sl_bar_slope, "ASD", "Tone")
asd_syllable <- stat_desc(iascl_sl_bar_slope, "ASD", "Syllable")
 
extract_line_graph_stat <- function (x, group, task) {
y <- x[c("nbr.val","mean", "SE.mean"),"value"]
dim(y) <- c(3, 1)
y <- t(y)
colnames(y) <- c("n", "mean", "std_error")
y <- data.frame(y)
y$group <- group
y$task <- task


return(y)
}

td_tone_stat <- extract_line_graph_stat(td_tone, "TD", "Tone")
td_syllable_stat <- extract_line_graph_stat(td_syllable, "TD", "Syllable")
asd_tone_stat <- extract_line_graph_stat(asd_tone, "ASD", "Tone")
asd_syllable_stat <- extract_line_graph_stat(asd_syllable, "ASD", "Syllable")

slope_stat <- rbind(td_tone_stat, td_syllable_stat, asd_tone_stat, asd_syllable_stat)

pd <- position_dodge(width = 0.2)

colnames(slope_stat)[colnames(slope_stat) == "group"] <- "Group"

library(plyr)
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:reshape':
## 
##     rename, round_any
slope_stat$task <-
revalue(slope_stat$task, c(#"entropy_children_lsl_entropy"="Letter",
                                     #"entropy_children_vsl_entropy"="Image",
                                     "Syllable"="Linguistic",
                                     "Tone"="Non-linguistic"))

pd <- position_dodge(width = 0.2)
slope_stat %>%
  ggplot(aes(x = task, y = mean, group = Group)) +
  geom_line(aes(linetype = Group, color = Group), position = pd, size = 1.8) +
   geom_errorbar(aes(ymin = mean - std_error, ymax = mean + std_error), width = .1, position = pd) +
     geom_point(aes(color = Group), size = 4, position = pd,show.legend = FALSE) + 
          geom_point(size = 3, color = "white", position = pd) +
      labs(title = "Behavioral Reaction Time Slope across Groups",
         x = "SL Task",  # Change x-axis label
         y = "Reaction Time Slope (arbitrary unit/trial)") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(
        plot.title = element_text(size=16, face="bold"),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text=element_text(size=12, face = "bold")
        ) +
  theme(legend.text=element_text(size=14, face="bold"),
        legend.title=element_text(size=15, face="bold")) +
   theme(
    panel.background = element_rect(fill = "white"),         # Set plot background to white
    legend.key  = element_rect(fill = "white"),              # Set legend item backgrounds to white
    axis.line.x = element_line(colour = "black", size = 1),  # Add line to x axis
    axis.line.y = element_line(colour = "black", size = 1)   # Add line to y axis
  ) +
  geom_text(aes(x= 1, y = 0.015, label = paste0("†")), size=7) 

# ggsave(
#    "behavioral_slope_across_group.png",
#     width = 15, height = 15, units = "cm"
# )

Behavioral Correlations

Heat Map

Heat Map for R, P, N for behavioral correlations

Within Task, Across Measures

#GGplot the heat map----------------------------------------------------------------------------------------------

ggplot(data = heat_map_corr_data, aes(x = column, y = row, fill = cor)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    midpoint = 0,
    limit = c(-1, 1),
    space = "Lab",
    name = "pearson\nCorrelation r"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 1,
    size = 12,
    hjust = 1
  )) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text (size = 12),
        axis.title.y = element_blank()) +
  coord_fixed() +
  geom_text(aes(label = cor), color = "black", size = 2)

ggplot(data = heat_map_corr_data, aes(x = column, y = row, fill = p)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    midpoint = 0,
    limit = c(0, 1),
    space = "Lab",
    name = "p-value"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 1,
    size = 12,
    hjust = 1
  )) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text (size = 12),
        axis.title.y = element_blank()) +
  coord_fixed() +
  geom_text(aes(label = p), color = "black", size = 2)

ggplot(data = heat_map_corr_data, aes(x = column, y = row, fill = n)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    midpoint = 0,
    limit = c(0, 33),
    space = "Lab",
    name = "n"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 1,
    size = 12,
    hjust = 1
  )) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text (size = 12),
        axis.title.y = element_blank()) +
  coord_fixed() +
  geom_text(aes(label = n), color = "black", size = 3.5)

Across Behavioral Tasks, Within Measure

ggplot(data = heat_map_corr_data_task, aes(x = column, y = row, fill = cor)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    midpoint = 0,
    limit = c(-1, 1),
    space = "Lab",
    name = "pearson\nCorrelation r"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 1,
    size = 12,
    hjust = 1
  )) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text (size = 12),
        axis.title.y = element_blank()) +
  coord_fixed() +
  geom_text(aes(label = cor), color = "black", size = 2)

ggplot(data = heat_map_corr_data_task, aes(x = column, y = row, fill = p)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    midpoint = 0,
    limit = c(0, 1),
    space = "Lab",
    name = "p-value"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 1,
    size = 12,
    hjust = 1
  )) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text (size = 12),
        axis.title.y = element_blank()) +
  coord_fixed() +
  geom_text(aes(label = p), color = "black", size = 2)

Beh Corr Plots

ggplot(data = iascl_all_completed,
       aes(y= gjt_std_score_gjt_web, 
           x = scq_total,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2") +
   guides(linetype = guide_legend("Group")) +
      labs(title = "Grammaticality Judgement Task and Social Communication Questionnaire",
         y = "Grammaticality Judgement Standardized Score",  # Change x-axis label
         x = "Social Communication Questionnaire") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  geom_text(aes(x = 20, y = 120, label = paste("rho =", "-0.76**")),color = "red",size=4)
## Warning: Removed 21 rows containing missing values (geom_point).

ggplot(data = iascl_all_completed,
       aes(y= a_prime, 
           x = scq_total,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2") +
   guides(linetype = guide_legend("Group")) +
      labs(title = "Grammaticality Judgement Task and Social Communication Questionnaire",
         y = "Grammaticality Judgement Raw Score",  # Change x-axis label
         x = "Social Communication Questionnaire") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  geom_text(aes(x = 20, y = 1, label = paste("rho =", "-0.65*")),color = "red",size=4)
## Warning: Removed 21 rows containing missing values (geom_point).

ggplot(data = iascl_all_completed,
       aes(x= a_prime, 
           y = scq_total,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2") +
   guides(linetype = guide_legend("Group")) +
      labs(title = "Grammaticality Judgement Task and Social Communication Questionnaire",
         x = "Grammaticality Judgement Raw Score",  # Change x-axis label
         y = "Social Communication Questionnaire") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  geom_text(aes(y = 20, x = 0.95, label = paste("rho =", "-0.65*")),color = "red",size=4)
## Warning: Removed 21 rows containing missing values (geom_point).

ggplot(data = iascl_all_completed,
       aes(x= a_prime, 
           y = scq_total,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2") +
   guides(linetype = guide_legend("Group")) +
      labs(title = "Grammaticality Judgement Task and Social Communication Questionnaire",
         x = "Grammaticality Judgement Raw Score",  # Change x-axis label
         y = "Social Communication Questionnaire") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  geom_text(aes(y = 20, x = 0.95, label = paste("rho =", "-0.65*")),color = "red",size=4)
## Warning: Removed 21 rows containing missing values (geom_point).

Ignore: ASD Beh Corr Plot

#MRI T-test between group

iascl_beh_gray_matter_td[,c(2:10, 13:23)] <- 
  lapply(iascl_beh_gray_matter_td[,c(2:10, 13:23)], function(x) {
  if (is.factor(x))
    as.numeric(as.character(x))
  else
    x
})


iascl_beh_gray_matter_asd[,c(2:10, 13:23)] <- 
  lapply(iascl_beh_gray_matter_asd[,c(2:10, 13:23)], function(x) {
  if (is.factor(x))
    as.numeric(as.character(x))
  else
    x
})


iascl_beh_cor_thick_td[,c(2:10, 13:23)] <- 
  lapply(iascl_beh_cor_thick_td[,c(2:10, 13:23)], function(x) {
  if (is.factor(x))
    as.numeric(as.character(x))
  else
    x
})


iascl_beh_cor_thick_asd[,c(2:10, 13:23)] <- 
  lapply(iascl_beh_cor_thick_asd[,c(2:10, 13:23)], function(x) {
  if (is.factor(x))
    as.numeric(as.character(x))
  else
    x
})


iascl_beh_gray_matter[,c(2:10, 13:23)] <- 
  lapply(iascl_beh_gray_matter[,c(2:10, 13:23)], function(x) {
  if (is.factor(x))
    as.numeric(as.character(x))
  else
    x
})


iascl_beh_cor_thick[,c(2:10, 13:23)] <- 
  lapply(iascl_beh_cor_thick[,c(2:10, 13:23)], function(x) {
  if (is.factor(x))
    as.numeric(as.character(x))
  else
    x
})


#T test between group
t_test_multi_pair <- 
  function(x,y){
  test <- t.test(x,y)

  data.frame(p_value = test$p.value,
             df = test$parameter,
             t_stat = test$statistic)
  }

# GM T-test between group
sapply(intersect(colnames(iascl_beh_gray_matter_td),colnames(iascl_beh_gray_matter_asd))[2:10], 
                 function(x) t_test_multi_pair(iascl_beh_gray_matter_td[,x], iascl_beh_gray_matter_asd[,x]))
##         LeftAngG   LeftIFG   LeftIFGorb LeftMFG   LeftMidAntTemp
## p_value 0.3709168  0.2490072 0.3570282  0.2734187 0.4028914     
## df      5.571482   5.095522  4.784659   3.941682  3.900906      
## t_stat  -0.9729214 -1.300917 -1.018772  -1.271432 -0.9373715    
##         LeftMidPostTemp LeftPostTemp RightMidAntTemp RightMidPostTemp
## p_value 0.2847049       0.6574654    0.3053455       0.1005669       
## df      6.152649        8.325675     3.941514        9.83297         
## t_stat  -1.171644       -0.4597477   -1.177003       -1.81211
# CT T-test between group
sapply(intersect(colnames(iascl_beh_cor_thick_td),colnames(iascl_beh_cor_thick_asd))[2:10], 
                 function(x) t_test_multi_pair(iascl_beh_cor_thick_td[,x], iascl_beh_cor_thick_asd[,x]))
##         LeftAngG  LeftIFG   LeftIFGorb LeftMFG    LeftMidAntTemp
## p_value 0.1295935 0.0905683 0.4910804  0.07046103 0.9654963     
## df      8.541535  10.5116   3.995009   4.960533   12.9725       
## t_stat  -1.677366 -1.863315 -0.7573    -2.29689   0.0441        
##         LeftMidPostTemp LeftPostTemp RightMidAntTemp RightMidPostTemp
## p_value 0.1718029       0.1839262    0.6907809       0.2909642       
## df      5.926682        7.287599     5.010944        7.925119        
## t_stat  -1.553923       -1.467877    -0.4216162      -1.131385
# ETIV T-test between group

No significant difference between TD and ASD for Cortical thickness or Gray Matter.

ETIV-Gray Matter Corr (Only 13 have ETIV + GM)

etiv <- 
read.csv("/Users/jojohu/Documents/Qlab/iascl/aseg_stats_etiv_measure.csv")

library("stringr")

etiv$Measure.volume <- str_extract(etiv$Measure.volume, "(?<=sub-)\\S{9}(?=/)")

etiv$Measure.volume 
##  [1] "blasta001" "blasta002" "blasta004" "blasta005" "blasta006"
##  [6] "blasta007" "blasta008" "blasta010" "blasta011" "blasta012"
## [11] "blasta013" "blasta015" "blasta017" "blasta018" "blasta020"
## [16] "blasta021" "blasta022" "blasta023" "blasta024" "blasta025"
## [21] "blasta026" "blasta027" "blasta028" "blasta029" "blastc009"
## [26] "blastc019" "blastc059" "blastc061" "blastc071" "blastc095"
## [31] "blastc126" "blastc147" "blastc165" "blastc168" "blastc206"
## [36] "blastc210" "blastc213" "blastc224"
partid1 <- substr(etiv$Measure.volume, start = 1, stop = 5)

partid2 <- substr(etiv$Measure.volume, start = 6, stop = 6)

partid3 <- substr(etiv$Measure.volume, start = 7, stop = 9)

etiv$part_id <- paste(partid1, partid2, partid3, sep = "_")

mri_etiv_gm_wide <- 
merge(mri_gray_matter_wide, 
      etiv[,c("part_id", "EstimatedTotalIntraCranialVol")],
      by.x = "part_id", by.y = "part_id",
      all.x = TRUE)

etiv_with_group <-
merge(iascl_beh_cor_thick, etiv[,c("part_id", "EstimatedTotalIntraCranialVol")],
      by.x = "part_id", by.y="part_id",
      all.x = TRUE)

#T test for ETIV
t.test(etiv_with_group[which(etiv_with_group$group == "TD"), "EstimatedTotalIntraCranialVol"], 
       etiv_with_group[which(etiv_with_group$group == "ASD"), "EstimatedTotalIntraCranialVol"])
## 
##  Welch Two Sample t-test
## 
## data:  etiv_with_group[which(etiv_with_group$group == "TD"), "EstimatedTotalIntraCranialVol"] and etiv_with_group[which(etiv_with_group$group == "ASD"), "EstimatedTotalIntraCranialVol"]
## t = -1.5014, df = 10.338, p-value = 0.1632
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -248036.55   47803.49
## sample estimates:
## mean of x mean of y 
##   1549337   1649453
mri_etiv_gm_corr <- rcorr(as.matrix(
  mri_etiv_gm_wide[,c(2:11)]), type = c("pearson"))


mri_etiv_gm_pearson_corr <- 
  flattenCorrMatrix(mri_etiv_gm_corr$r, mri_etiv_gm_corr$P, mri_etiv_gm_corr$n)


mri_etiv_gm_pearson_corr[,c(3:5)] <- 
  round(mri_etiv_gm_pearson_corr[,c(3:5)], 3)

mri_etiv_gm_pearson_corr <-
  mri_etiv_gm_pearson_corr[which(mri_etiv_gm_pearson_corr$p <= 0.07),]


print(mri_etiv_gm_pearson_corr)
##                 row                        column   cor     p  n
## 1          LeftAngG                       LeftIFG 0.813 0.000 16
## 2          LeftAngG                    LeftIFGorb 0.849 0.000 16
## 3           LeftIFG                    LeftIFGorb 0.847 0.000 16
## 4          LeftAngG                       LeftMFG 0.688 0.003 16
## 5           LeftIFG                       LeftMFG 0.858 0.000 16
## 6        LeftIFGorb                       LeftMFG 0.695 0.003 16
## 8           LeftIFG                LeftMidAntTemp 0.464 0.070 16
## 9        LeftIFGorb                LeftMidAntTemp 0.559 0.024 16
## 10          LeftMFG                LeftMidAntTemp 0.714 0.002 16
## 11         LeftAngG               LeftMidPostTemp 0.706 0.002 16
## 12          LeftIFG               LeftMidPostTemp 0.733 0.001 16
## 13       LeftIFGorb               LeftMidPostTemp 0.639 0.008 16
## 14          LeftMFG               LeftMidPostTemp 0.577 0.019 16
## 16         LeftAngG                  LeftPostTemp 0.643 0.007 16
## 17          LeftIFG                  LeftPostTemp 0.833 0.000 16
## 18       LeftIFGorb                  LeftPostTemp 0.673 0.004 16
## 19          LeftMFG                  LeftPostTemp 0.706 0.002 16
## 21  LeftMidPostTemp                  LeftPostTemp 0.655 0.006 16
## 25          LeftMFG               RightMidAntTemp 0.480 0.060 16
## 26   LeftMidAntTemp               RightMidAntTemp 0.619 0.011 16
## 28     LeftPostTemp               RightMidAntTemp 0.530 0.035 16
## 29         LeftAngG              RightMidPostTemp 0.515 0.041 16
## 30          LeftIFG              RightMidPostTemp 0.631 0.009 16
## 32          LeftMFG              RightMidPostTemp 0.655 0.006 16
## 34  LeftMidPostTemp              RightMidPostTemp 0.483 0.058 16
## 35     LeftPostTemp              RightMidPostTemp 0.613 0.012 16
## 36  RightMidAntTemp              RightMidPostTemp 0.496 0.051 16
## 37         LeftAngG EstimatedTotalIntraCranialVol 0.696 0.008 13
## 38          LeftIFG EstimatedTotalIntraCranialVol 0.694 0.009 13
## 39       LeftIFGorb EstimatedTotalIntraCranialVol 0.773 0.002 13
## 40          LeftMFG EstimatedTotalIntraCranialVol 0.568 0.043 13
## 41   LeftMidAntTemp EstimatedTotalIntraCranialVol 0.555 0.049 13
## 42  LeftMidPostTemp EstimatedTotalIntraCranialVol 0.581 0.037 13
## 45 RightMidPostTemp EstimatedTotalIntraCranialVol 0.700 0.008 13

ETIV not significantly different between group.

LeftAngG, LeftIFG, LeftIFGorb, LeftMFG, LeftMidAntTemp are correlated with ETIV.

ETIV-Cortical Thickness Correlation (only 13 have ETIV plus Cortical Thickness)

mri_etiv_ct_wide <- 
merge(average_cor_thick_wide, 
      etiv[,c("part_id", "EstimatedTotalIntraCranialVol")],
      by.x = "part_id", by.y = "part_id",
      all.x = TRUE)


mri_etiv_ct_corr <- rcorr(as.matrix(
  mri_etiv_ct_wide[,c(2:11)]), type = c("pearson"))


mri_etiv_ct_pearson_corr <- 
  flattenCorrMatrix(mri_etiv_ct_corr$r, mri_etiv_ct_corr$P, mri_etiv_ct_corr$n)


mri_etiv_ct_pearson_corr[,c(3:5)] <- 
  round(mri_etiv_ct_pearson_corr[,c(3:5)], 3)

mri_etiv_ct_pearson_corr <-
  mri_etiv_ct_pearson_corr[which(mri_etiv_ct_pearson_corr$p <= 0.07),]


print(mri_etiv_ct_pearson_corr)
##                row           column   cor     p  n
## 1         LeftAngG          LeftIFG 0.499 0.049 16
## 3          LeftIFG       LeftIFGorb 0.659 0.006 16
## 4         LeftAngG          LeftMFG 0.635 0.008 16
## 5          LeftIFG          LeftMFG 0.772 0.000 16
## 12         LeftIFG  LeftMidPostTemp 0.754 0.001 16
## 13      LeftIFGorb  LeftMidPostTemp 0.469 0.067 16
## 14         LeftMFG  LeftMidPostTemp 0.683 0.004 16
## 15  LeftMidAntTemp  LeftMidPostTemp 0.581 0.018 16
## 16        LeftAngG     LeftPostTemp 0.552 0.027 16
## 17         LeftIFG     LeftPostTemp 0.587 0.017 16
## 19         LeftMFG     LeftPostTemp 0.733 0.001 16
## 21 LeftMidPostTemp     LeftPostTemp 0.488 0.055 16
## 29        LeftAngG RightMidPostTemp 0.518 0.040 16
## 30         LeftIFG RightMidPostTemp 0.522 0.038 16
## 31      LeftIFGorb RightMidPostTemp 0.532 0.034 16

Cortical thickness is not correlated with ETIV.

MRI-Behavioral Corr (ASD + TD without behavioral outliers)

Gray Matter-Beh Corr (ASD + TD)

iascl_beh_gray_matter_all <- rcorr(as.matrix(
  iascl_beh_gray_matter[,c(2:10, 13:23)]), type = c("pearson"))


iascl_beh_gray_matter_pearson_all <- 
  flattenCorrMatrix(iascl_beh_gray_matter_all$r, iascl_beh_gray_matter_all$P, iascl_beh_gray_matter_all$n)


iascl_beh_gray_matter_pearson_all[,c(3:5)] <- 
  round(iascl_beh_gray_matter_pearson_all[,c(3:5)], 3)

iascl_beh_gray_matter_pearson_all <-
  iascl_beh_gray_matter_pearson_all[which(iascl_beh_gray_matter_pearson_all$p <= 0.05),]


print(iascl_beh_gray_matter_pearson_all)
##                                  row                             column
## 1                           LeftAngG                            LeftIFG
## 2                           LeftAngG                         LeftIFGorb
## 3                            LeftIFG                         LeftIFGorb
## 4                           LeftAngG                            LeftMFG
## 5                            LeftIFG                            LeftMFG
## 6                         LeftIFGorb                            LeftMFG
## 9                         LeftIFGorb                     LeftMidAntTemp
## 10                           LeftMFG                     LeftMidAntTemp
## 11                          LeftAngG                    LeftMidPostTemp
## 12                           LeftIFG                    LeftMidPostTemp
## 13                        LeftIFGorb                    LeftMidPostTemp
## 14                           LeftMFG                    LeftMidPostTemp
## 16                          LeftAngG                       LeftPostTemp
## 17                           LeftIFG                       LeftPostTemp
## 18                        LeftIFGorb                       LeftPostTemp
## 19                           LeftMFG                       LeftPostTemp
## 21                   LeftMidPostTemp                       LeftPostTemp
## 26                    LeftMidAntTemp                    RightMidAntTemp
## 28                      LeftPostTemp                    RightMidAntTemp
## 29                          LeftAngG                   RightMidPostTemp
## 30                           LeftIFG                   RightMidPostTemp
## 32                           LeftMFG                   RightMidPostTemp
## 35                      LeftPostTemp                   RightMidPostTemp
## 41                    LeftMidAntTemp              gjt_std_score_gjt_web
## 44                   RightMidAntTemp              gjt_std_score_gjt_web
## 45                  RightMidPostTemp              gjt_std_score_gjt_web
## 91  accuracy_children_tsl_accuracies       entropy_children_ssl_entropy
## 104 accuracy_children_tsl_accuracies       entropy_children_tsl_entropy
## 106                         LeftAngG          rt_children_ssl_indiv_rts
## 108                       LeftIFGorb          rt_children_ssl_indiv_rts
## 109                          LeftMFG          rt_children_ssl_indiv_rts
## 111                  LeftMidPostTemp          rt_children_ssl_indiv_rts
## 126                  LeftMidPostTemp          rt_children_tsl_indiv_rts
## 142                  LeftMidPostTemp slope_children_ssl_indiv_rts_slope
## 153        rt_children_tsl_indiv_rts slope_children_ssl_indiv_rts_slope
## 154                         LeftAngG slope_children_tsl_indiv_rts_slope
## 155                          LeftIFG slope_children_tsl_indiv_rts_slope
## 156                       LeftIFGorb slope_children_tsl_indiv_rts_slope
## 157                          LeftMFG slope_children_tsl_indiv_rts_slope
## 159                  LeftMidPostTemp slope_children_tsl_indiv_rts_slope
## 160                     LeftPostTemp slope_children_tsl_indiv_rts_slope
## 169        rt_children_ssl_indiv_rts slope_children_tsl_indiv_rts_slope
## 170        rt_children_tsl_indiv_rts slope_children_tsl_indiv_rts_slope
## 172                         LeftAngG                  age_at_neuro_year
## 185     entropy_children_ssl_entropy                  age_at_neuro_year
##        cor     p  n
## 1    0.813 0.000 16
## 2    0.849 0.000 16
## 3    0.847 0.000 16
## 4    0.688 0.003 16
## 5    0.858 0.000 16
## 6    0.695 0.003 16
## 9    0.559 0.024 16
## 10   0.714 0.002 16
## 11   0.706 0.002 16
## 12   0.733 0.001 16
## 13   0.639 0.008 16
## 14   0.577 0.019 16
## 16   0.643 0.007 16
## 17   0.833 0.000 16
## 18   0.673 0.004 16
## 19   0.706 0.002 16
## 21   0.655 0.006 16
## 26   0.619 0.011 16
## 28   0.530 0.035 16
## 29   0.515 0.041 16
## 30   0.631 0.009 16
## 32   0.655 0.006 16
## 35   0.613 0.012 16
## 41  -0.938 0.018  5
## 44  -0.917 0.028  5
## 45  -0.964 0.008  5
## 91  -0.732 0.010 11
## 104 -0.730 0.007 12
## 106 -0.597 0.031 13
## 108 -0.607 0.028 13
## 109 -0.573 0.041 13
## 111 -0.718 0.006 13
## 126  0.596 0.041 12
## 142 -0.602 0.030 13
## 153 -0.612 0.045 11
## 154  0.739 0.006 12
## 155  0.739 0.006 12
## 156  0.604 0.037 12
## 157  0.665 0.018 12
## 159  0.876 0.000 12
## 160  0.795 0.002 12
## 169 -0.619 0.042 11
## 170  0.659 0.020 12
## 172 -0.524 0.037 16
## 185  0.732 0.004 13
# write.csv(iascl_beh_gray_matter_pearson_all, "iascl_beh_gray_matter_pearson_all.csv")

Ignore: Gray Matter-Beh Corr TD Only

Partial Correlations Controling for Age and IQ

partial_corr_data <- 
merge(iascl_beh_gray_matter, drive_checklist[,c("Participant.ID", "kbit_matrices_std")],
      by.x = "part_id", by.y = "Participant.ID",
      all.x = TRUE)

partial_corr_data$kbit_matrices_std <- as.numeric(as.character(partial_corr_data$kbit_matrices_std))

partial_corr_data_ssl <- subset(partial_corr_data, !is.na(rt_children_ssl_indiv_rts)) 
partial_corr_data_tsl <- subset(partial_corr_data, !is.na(rt_children_tsl_indiv_rts)) 

library("ppcor")
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
partial_corr_tsl_result <- 
  pcor.test(partial_corr_data_tsl$LeftMidPostTemp,
            partial_corr_data_tsl$rt_children_tsl_indiv_rts,
            partial_corr_data_tsl[,c("age_at_neuro_year","kbit_matrices_std")],method="pearson")
partial_corr_tsl_result <- describe(partial_corr_tsl_result)

partial_corr_ssl_result <- 
  pcor.test(partial_corr_data_ssl$LeftMidPostTemp,
            partial_corr_data_ssl$rt_children_ssl_indiv_rts,
            partial_corr_data_ssl[,c("age_at_neuro_year","kbit_matrices_std")],method="pearson")
partial_corr_ssl_result <- describe(partial_corr_ssl_result)






library("psych")
## 
## Attaching package: 'psych'
## The following object is masked from 'package:Hmisc':
## 
##     describe
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
partial_cor_r <- partial.r(data=partial_corr_data[,-c(1,11, 12)], 
                           x=c("LeftMidPostTemp", "rt_children_ssl_indiv_rts"),
                           y=c("age_at_neuro_year", "kbit_matrices_std"), method = "pearson", use="pairwise")
 
cp <- corr.p(partial_cor_r,n=11)
print(cp, short = F)
## Call:corr.p(r = partial_cor_r, n = 11)
## Correlation matrix 
## partial correlations 
##                           LeftMidPostTemp rt_children_ssl_indiv_rts
## LeftMidPostTemp                      1.00                     -0.71
## rt_children_ssl_indiv_rts           -0.71                      1.00
## Sample Size 
## [1] 11
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
## partial correlations 
##                           LeftMidPostTemp rt_children_ssl_indiv_rts
## LeftMidPostTemp                      0.00                      0.01
## rt_children_ssl_indiv_rts            0.01                      0.00
## 
##  Confidence intervals based upon normal theory.  To get bootstrapped values, try cor.ci
##             lower     r upper    p
## LfMPT-r____ -0.92 -0.71  -0.2 0.01
partial_cor_r_tsl <- partial.r(data=partial_corr_data[,-c(1,11, 12)], 
                           x=c("LeftMidPostTemp", "rt_children_tsl_indiv_rts"),
                           y=c("age_at_neuro_year", "kbit_matrices_std"), method = "pearson", use="pairwise")
 
cp_tsl <- corr.p(partial_cor_r_tsl,n=10)
print(cp_tsl, short = F)
## Call:corr.p(r = partial_cor_r_tsl, n = 10)
## Correlation matrix 
## partial correlations 
##                           LeftMidPostTemp rt_children_tsl_indiv_rts
## LeftMidPostTemp                      1.00                      0.74
## rt_children_tsl_indiv_rts            0.74                      1.00
## Sample Size 
## [1] 10
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
## partial correlations 
##                           LeftMidPostTemp rt_children_tsl_indiv_rts
## LeftMidPostTemp                      0.00                      0.01
## rt_children_tsl_indiv_rts            0.01                      0.00
## 
##  Confidence intervals based upon normal theory.  To get bootstrapped values, try cor.ci
##             lower    r upper    p
## LfMPT-r____  0.21 0.74  0.93 0.01

Ignore: Gray Matter-Beh Corr ASD Only

Avg Cortical Thickness-Beh Corr (ASD + TD)

iascl_beh_cor_thick_pearson_all <- rcorr(as.matrix(
  iascl_beh_cor_thick[,c(2:10, 13:23)]), type = c("pearson"))


iascl_beh_cor_thick_pearson_all <- 
  flattenCorrMatrix(iascl_beh_cor_thick_pearson_all$r,
                    iascl_beh_cor_thick_pearson_all$P, 
                    iascl_beh_cor_thick_pearson_all$n)


iascl_beh_cor_thick_pearson_all[,c(3:5)] <- 
  round(iascl_beh_cor_thick_pearson_all[,c(3:5)], 3)

iascl_beh_cor_thick_pearson_all <-
  iascl_beh_cor_thick_pearson_all[which(iascl_beh_cor_thick_pearson_all$p <= 0.05),]


print(iascl_beh_cor_thick_pearson_all)
##                                  row                             column
## 1                           LeftAngG                            LeftIFG
## 3                            LeftIFG                         LeftIFGorb
## 4                           LeftAngG                            LeftMFG
## 5                            LeftIFG                            LeftMFG
## 12                           LeftIFG                    LeftMidPostTemp
## 14                           LeftMFG                    LeftMidPostTemp
## 15                    LeftMidAntTemp                    LeftMidPostTemp
## 16                          LeftAngG                       LeftPostTemp
## 17                           LeftIFG                       LeftPostTemp
## 19                           LeftMFG                       LeftPostTemp
## 29                          LeftAngG                   RightMidPostTemp
## 30                           LeftIFG                   RightMidPostTemp
## 31                        LeftIFGorb                   RightMidPostTemp
## 91  accuracy_children_tsl_accuracies       entropy_children_ssl_entropy
## 104 accuracy_children_tsl_accuracies       entropy_children_tsl_entropy
## 111                  LeftMidPostTemp          rt_children_ssl_indiv_rts
## 153        rt_children_tsl_indiv_rts slope_children_ssl_indiv_rts_slope
## 169        rt_children_ssl_indiv_rts slope_children_tsl_indiv_rts_slope
## 170        rt_children_tsl_indiv_rts slope_children_tsl_indiv_rts_slope
## 185     entropy_children_ssl_entropy                  age_at_neuro_year
##        cor     p  n
## 1    0.499 0.049 16
## 3    0.659 0.006 16
## 4    0.635 0.008 16
## 5    0.772 0.000 16
## 12   0.754 0.001 16
## 14   0.683 0.004 16
## 15   0.581 0.018 16
## 16   0.552 0.027 16
## 17   0.587 0.017 16
## 19   0.733 0.001 16
## 29   0.518 0.040 16
## 30   0.522 0.038 16
## 31   0.532 0.034 16
## 91  -0.732 0.010 11
## 104 -0.730 0.007 12
## 111 -0.633 0.020 13
## 153 -0.612 0.045 11
## 169 -0.619 0.042 11
## 170  0.659 0.020 12
## 185  0.732 0.004 13
# write.csv(iascl_beh_cor_thick_pearson_all, "iascl_beh_cor_thick_pearson_all.csv")

Ignore: Avg Cortical Thickness-Beh Corr TD Only

Ignore: Avg Cortical Thickness-Beh Corr ASD Only

MRI Gray Matter Corr Plots

SSL Mean RT: Negative Corr between SSL Mean RT and Gray Matter Volume

ggplot(data = subset(iascl_beh_gray_matter, !is.na(rt_children_ssl_indiv_rts)),
       aes(x=rt_children_ssl_indiv_rts, 
           y = LeftAngG,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

ggplot(data = subset(iascl_beh_gray_matter, !is.na(rt_children_ssl_indiv_rts)),
       aes(x=rt_children_ssl_indiv_rts, 
           y = LeftMidPostTemp,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

iascl_beh_mri_corr_plot <-
iascl_beh_gray_matter[, c("part_id",
                          "LeftMidPostTemp", 
                          "slope_children_ssl_indiv_rts_slope",
                          "slope_children_tsl_indiv_rts_slope",
                          "group")]


iascl_beh_mri_corr_plot <- melt(iascl_beh_mri_corr_plot, id.vars=c("part_id", "group", "LeftMidPostTemp"))

iascl_beh_mri_corr_plot[which(iascl_beh_mri_corr_plot$part_id == "blast_c_215"),
                        "group"] <- "TD"


iascl_beh_mri_corr_plot$variable <-
revalue(iascl_beh_mri_corr_plot$variable, c(#"accuracy_children_lsl_accuracies"="Letter",
                                     #"accuracy_children_vsl_accuracies"="Image",
                                     "slope_children_tsl_indiv_rts_slope"="Non-linguistic",
                                     "slope_children_ssl_indiv_rts_slope"="Linguistic"))

colnames(iascl_beh_mri_corr_plot)[colnames(iascl_beh_mri_corr_plot) == "variable"] <- "Task"
colnames(iascl_beh_mri_corr_plot)[colnames(iascl_beh_mri_corr_plot) == "group"] <- "Group"

iascl_beh_mri_corr_plot$legend <- paste(iascl_beh_mri_corr_plot$Group, iascl_beh_mri_corr_plot$Task)

library("stringr")
iascl_beh_mri_corr_plot$unique_id <- str_extract(iascl_beh_mri_corr_plot$part_id, "(?<=blast_c_)\\S+")

# write.csv(iascl_beh_mri_corr_plot, "iascl_beh_mri_corr_plot_data.csv")

library(ggpubr)
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:pastecs':
## 
##     extract
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
## 
##     mutate
library(plyr)

iascl_beh_mri_corr_plot %>%
ggplot(
       aes(x=value, 
           y = LeftMidPostTemp,
           color = legend,
           shape = legend
           )) +
  theme(legend.title = element_blank()) +
  geom_point(size = 4) +
   scale_shape_manual(values=c(16, 16, 17, 17))+
  scale_color_manual(values = c("#756bb1", "#fec44f", "#756bb1", "#fec44f")) +
   labs(title = "Grey Matter Volume Relationship with Behavioral Task \n",
         y = bquote(bold("Left Posterior Middle Temporal Lobe GMV"~(mm^3))),  # Change x-axis label
         x = "Reaction Time Slope (arbitrary unit/trial)") +
  theme(plot.title = element_text(hjust = 0.23)) +
  theme(
        plot.title = element_text(size=16, face="bold"),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text=element_text(size=12, face = "bold")
        ) +
  theme(legend.text=element_text(size=14, face="bold")) +
   theme(
    panel.background = element_rect(fill = "white"),         # Set plot background to white
    legend.key  = element_rect(fill = "white"),              # Set legend item backgrounds to white
    axis.line.x = element_line(colour = "black", size = 1),  # Add line to x axis
    axis.line.y = element_line(colour = "black", size = 1)   # Add line to y axis
  ) +
  geom_smooth(data = subset(iascl_beh_mri_corr_plot, Task == "Linguistic"),
              aes(x=value, y=LeftMidPostTemp),
              method=lm, se=FALSE, show.legend = F, inherit.aes = F, color = "#756bb1") +
    geom_smooth(data = subset(iascl_beh_mri_corr_plot, Task == "Non-linguistic"),
              aes(x=value, y=LeftMidPostTemp),
              method=lm, se=FALSE, show.legend = F, inherit.aes = F, color = "#fec44f") +
  stat_cor(data = subset(iascl_beh_mri_corr_plot, Task == "Linguistic"),
              aes(x=value, y=LeftMidPostTemp), method = "pearson", label.x = -0.05, label.y = 1590,
           inherit.aes = F, color = "#756bb1", size = 5) +
    stat_cor(data = subset(iascl_beh_mri_corr_plot, Task == "Non-linguistic"),
              aes(x=value, y=LeftMidPostTemp), method = "pearson", label.x = -0.05, label.y = 1550, 
           inherit.aes = F, color = "#fec44f", size = 5) #+
## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing non-finite values (stat_cor).
## Warning: Removed 4 rows containing non-finite values (stat_cor).
## Warning: Removed 7 rows containing missing values (geom_point).

   #geom_text(aes(label = unique_id),hjust=0.5, vjust=0.5, color = "black")

# ggsave("iascl_beh_mri_slope_corr_plot.png",
#        width = 15, height = 15, units = "cm")

MRI Behavioral correlation plot with text labels

library(ggpubr)
library(plyr)

iascl_beh_mri_corr_plot_text <- iascl_beh_mri_corr_plot

iascl_beh_mri_corr_plot_text$legend <- 
  factor(iascl_beh_mri_corr_plot_text$legend,levels = c("ASD Linguistics", 
                                              "TD Linguistics", 
                                              "ASD Nonlinguistics", 
                                              "TD Nonlinguistics"))

iascl_sl_bar_acc <-
  iascl_sl_bar_acc[order(iascl_sl_bar_acc$variable), ]



iascl_beh_mri_corr_plot_text %>%
ggplot(
       aes(x=value, 
           y = LeftMidPostTemp,
           color = legend
           #shape = legend
           )) +
  theme(legend.title = element_blank()) +
  scale_fill_discrete(labels = c("A", "B", "Linguistics", "Non-Linguistics")) +
  geom_point(size = 4) +
   scale_shape_manual(values=c(16, 16))+
  #scale_color_manual(values = c("#756bb1", "#fec44f", "#756bb1", "#fec44f")) +
  #scale_color_manual(values = c("#756bb1", "#fec44f")) +
   scale_colour_manual(labels = c("ASD", "ASD", "TD", "TD"),
     values = c("#756bb1", "#756bb1", "#fec44f", "#fec44f")) +
    #                   guide = guide_legend(override.aes = list(
     #                    linetype = c(rep("blank", 2), "solid", "dotted"), 
      #                   shape = c(rep(16, 2), NA, NA)))) +
   labs(title = "Grey Matter Volume Relationship with Behavioral Task \n",
         y = bquote(bold("Left Posterior Middle Temporal Lobe GMV"~(mm^3))),  # Change x-axis label
         x = "Reaction Time Slope (arbitrary unit/trial)") +
  theme(plot.title = element_text(hjust = 0.23)) +
  theme(
        plot.title = element_text(size=16, face="bold"),
        axis.title.x = element_text(size=14, face="bold"),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text=element_text(size=12, face = "bold")
        ) +
  theme(legend.text=element_text(size=14, face="bold")) +
   theme(
    panel.background = element_rect(fill = "white"),         # Set plot background to white
    legend.key  = element_rect(fill = "white"),              # Set legend item backgrounds to white
    axis.line.x = element_line(colour = "black", size = 1),  # Add line to x axis
    axis.line.y = element_line(colour = "black", size = 1)   # Add line to y axis
  ) +
  geom_smooth(data = subset(iascl_beh_mri_corr_plot, Task == "Linguistic"),
              aes(x=value, y=LeftMidPostTemp), linetype = "solid",
              method=lm, se=FALSE, show.legend = T, inherit.aes = F, color = "black") +
    geom_smooth(data = subset(iascl_beh_mri_corr_plot, Task == "Non-linguistic"),
              aes(x=value, y=LeftMidPostTemp), linetype = "dotted",
              method=lm, se=FALSE, show.legend = T, inherit.aes = F, color = "black") +
  #stat_cor(data = subset(iascl_beh_mri_corr_plot, Task == "Linguistic"),
  #            aes(x=value, y=LeftMidPostTemp), method = "pearson", label.x = -0.05, label.y = 1590,
   #        inherit.aes = F, color = "#756bb1", size = 5) +
    #stat_cor(data = subset(iascl_beh_mri_corr_plot, Task == "Non-linguistic"),
    #          aes(x=value, y=LeftMidPostTemp), method = "pearson", label.x = -0.05, label.y = 1550, 
     #      inherit.aes = F, color = "#fec44f", size = 5) +
   geom_text(aes(label = unique_id),hjust=0.8, vjust=0.8, color = "black", face = "bold", size =3) +
  geom_text(aes(y = 1590, x = -0.05, label = paste("R =", "0.74*"), fontface = 3),color = "black",size=4) +
  geom_text(aes(y = 1550, x = -0.05, label = paste("R =", "-0.71*"),  fontface = 3),color = "black",size=4)
## Warning: Ignoring unknown parameters: face
## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 32 rows containing missing values (geom_point).
## Warning: Removed 7 rows containing missing values (geom_text).

# ggsave("iascl_beh_mri_slope_corr_plot_with_text.png",
#        width = 15, height = 15, units = "cm")

TSL Mean RT: Positive Corr between TSL Mean RT and Gray Matter Volume

ggplot(data = subset(iascl_beh_gray_matter, !is.na(rt_children_tsl_indiv_rts)),
       aes(x=rt_children_tsl_indiv_rts, 
           y = LeftMidPostTemp,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

TSL Slope: Positive Corr between TSL RT Slope and Gray Matter Volume

ggplot(data = subset(iascl_beh_gray_matter, !is.na(slope_children_tsl_indiv_rts_slope)),
       aes(x=slope_children_tsl_indiv_rts_slope, 
           y = LeftAngG,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

ggplot(data = subset(iascl_beh_gray_matter, !is.na(slope_children_tsl_indiv_rts_slope)),
       aes(x=slope_children_tsl_indiv_rts_slope, 
           y = LeftIFG,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

ggplot(data = subset(iascl_beh_gray_matter, !is.na(slope_children_tsl_indiv_rts_slope)),
       aes(x=slope_children_tsl_indiv_rts_slope, 
           y = LeftIFGorb,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

ggplot(data = subset(iascl_beh_gray_matter, !is.na(slope_children_tsl_indiv_rts_slope)),
       aes(x=slope_children_tsl_indiv_rts_slope, 
           y = LeftMFG,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

ggplot(data = subset(iascl_beh_gray_matter, !is.na(slope_children_tsl_indiv_rts_slope)),
       aes(x=slope_children_tsl_indiv_rts_slope, 
           y = LeftMidPostTemp,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

ggplot(data = subset(iascl_beh_gray_matter, !is.na(slope_children_tsl_indiv_rts_slope)),
       aes(x=slope_children_tsl_indiv_rts_slope, 
           y = LeftPostTemp,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

MRI Cortical Thickness Corr Plots

SSL Mean RT: Negative Corr between SSL Mean RT and Cortical Thickness

ggplot(data = subset(iascl_beh_cor_thick, !is.na(rt_children_ssl_indiv_rts)),
       aes(x=rt_children_ssl_indiv_rts, 
           y = LeftMidPostTemp,
           color = group)) +
  geom_point() + 
  scale_fill_brewer(palette = "Set2")

Gender MRI

drive_list <- read.csv("/Users/jojohu/Documents/Qlab/iascl/drive_checklist.csv")

drive_list <- drive_list[,c("Participant.ID", "Group", "Gender", "date.of.birth", 
            "kbit_matrices_raw", "kbit_matrices_std", "date_of_lab_assessment", "date_of_mri" )]

demo_mri <-
merge(iascl_beh_cor_thick, drive_list,
      by.x = "part_id", by.y = "Participant.ID",
      all.x = TRUE)


demo_mri_long <- melt(demo_mri, id.vars = c("part_id", "group"))
## Warning: attributes are not identical across measure variables; they will
## be dropped
demo_mri_gender_long <- demo_mri_long[which(demo_mri_long$variable == "Gender"),]

demo_mri_gender_long <- subset(demo_mri_gender_long, part_id != "blast_c_215")

add_row <- data.frame(matrix(ncol = 4, nrow = 1))
colnames(add_row) <- colnames(demo_mri_gender_long)
add_row$part_id <- "blast_c_215"
add_row$group <- "TD"
add_row$variable <- "Gender"
add_row$value <- "F"

demo_mri_gender_long <- rbind(demo_mri_gender_long, add_row)

gender_mri_wide <- cast(demo_mri_gender_long, group~value, length)


print(gender_mri_wide)
##   group F M
## 1   ASD 1 3
## 2    TD 9 3
print(chisq.test(gender_mri_wide))
## Warning in chisq.test(gender_mri_wide): Chi-squared approximation may be
## incorrect
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  gender_mri_wide
## X-squared = 1.4222, df = 1, p-value = 0.233

Gender is matched in the MRI pariticipants.

KBITT MRI

demo_mri_kbitt_long <- demo_mri_long[which(
  demo_mri_long$variable == "kbit_matrices_raw" | 
    demo_mri_long$variable == "kbit_matrices_std"),]

demo_mri_kbitt_long$value <- as.numeric(as.character(demo_mri_kbitt_long$value))

demo_mri_kbitt_long[which(demo_mri_kbitt_long$part_id == "blast_c_215"), "group"] <- "TD"

library("reshape")

kbitt_mri_group_mean <- cast(demo_mri_kbitt_long, group~variable, "mean")

print(kbitt_mri_group_mean)
##   group kbit_matrices_raw kbit_matrices_std
## 1   ASD          30.50000          112.7500
## 2    TD          31.58333          111.3333

To Do:

-demo info for MRI - Age for in lab MRI appending. -Ancova for gender match - Done -Ancova for Whole brain? Still? No corr between etiv and CT and only one region significant in GM. -make the t-test results into a nice format please